subroutine movenp(cc,Energy_in)
use type
use var
use parallel
implicit none
include 'mpif.h'
integer :: i=0, cc, ID_ac, k,j, label
real(dp), intent(in) :: Energy_in(5)
real(dp) :: EE, E_ac(5), all_E(size), oldE, dE(size), P(size), P2(size), PP, rand
!integer, dimension(:,:,:), allocatable :: temp_po
logical :: flag, accept

if (cc/=-1) then
if (cc>1) then
   oldE = Energy_sycle(cc-1)
   
   EE = Energy_in(5)
   CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
   CALL  MPI_ALLGATHER(EE,1,MPI_DOUBLE_PRECISION,all_E,1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,ierr)
   dE(:) = all_E(:) - oldE
   dE = -beta * dE * Escale / kbt

   ! Max(1,e^dE)
   Do i=1,size
      if (dE(i)>0) dE(i)=0
   enddo
   P  = exp(dE)
   P  = P/sum(P)
   
   P2(1) = P(1)
   Do i=2,size
      P2(i) = P2(i-1)+P(i)
   enddo

   call random_number(rand)
   PP = 0
   Do i=1,size
      if (PP<=rand.and.rand<P2(i)) ID_ac = i - 1
      PP = P2(i)
   enddo
   
   E_ac = Energy_in
   CALL MPI_BCAST(ID_ac, 1, MPI_INTEGER, master, MPI_COMM_WORLD, IERR)
   CALL MPI_BCAST(NP_position, 3*number_NP, MPI_REAL, ID_ac, MPI_COMM_WORLD, IERR)
   CALL MPI_BCAST(E_ac, 5, MPI_DOUBLE_PRECISION, ID_ac, MPI_COMM_WORLD, IERR)
   CALL MPI_BCAST(q, 6*tnodes, MPI_DOUBLE_PRECISION, ID_ac, MPI_COMM_WORLD, IERR) !
   
   NP_position_back = NP_position ! backup for traj
   trj_E(1:5,1) = E_ac ! Save Energy for traj file
   call wtraj(cc,0,0)
      
   Energy_sycle(cc) = E_ac(5)
else
   Energy_sycle(cc) = Energy_in(5)
   ID_ac = 0
   E_ac = Energy_in
   CALL MPI_BCAST(q, 6*tnodes, MPI_DOUBLE_PRECISION, ID_ac, MPI_COMM_WORLD, IERR) !
   
   NP_position_back = NP_position ! backup for traj
   trj_E(1:5,1) = E_ac ! Save Energy for traj file
   call wtraj(cc,0,0)
endif


if (ID == master) then
   write(*,*) '############################################################'
   write(*,*) "# In cycle",cc,"the motion of",ID_ac+1,"was accepted #"
   write(*,'(5(f15.6))')  E_ac
   write(*,*)'############################################################'
   write(222,*) '############################################################'
   write(222,*) "# In cycle",cc,"the motion of",ID_ac+1,"was accepted #"
   write(222,'(5(f15.6))')  E_ac
   write(222,*)'############################################################'
   
   write(900,'(I5,2X,4(f15.6),2(f20.6))') cc, E_ac,Escale*E_ac(5)/kbT
   write(999,*)number_NP
   write(999,*)
   Do i=1,number_NP
      write(999,*)"O",NP_position(:,i)
   enddo
   open(77,file='restart')
   write(77,*)number_NP
   write(77,*)
   Do i=1,number_NP
      write(77,*)"O",NP_position(:,i)
   enddo
   close(77)
endif

endif
if (sycle>1.and.cc<sycle) call move
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)

call initial(cc)
call make_np()
  
  bulk = .false.
  nodes = 0
  do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
           if (number_NP>0) then
           if (channel) then
              if (boundary(i,j,k)==.false..and.boundary_NP(i,j,k)==.false..and.NP(i,j,k)==.false.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           elseif (droplet) then
              if (boundary(i,j,k)==.false..and.boundary_NP(i,j,k)==.false..and.NP(i,j,k)==.false..and.drop(i,j,k)==.true.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           endif
           else
           if (channel) then
              if (boundary(i,j,k)==.false.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           elseif (droplet) then
              if (boundary(i,j,k)==.false..and.drop(i,j,k)==.true.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           endif
           endif
        enddo
     enddo
  enddo

! Calculate q6
q(6,:) = -q(1,:) - q(4,:)

contains
subroutine move()
use var
use type
use parallel
implicit none
integer :: R, NP, i, delta, x, y
real :: pos(3)
real(dp) :: rand, rand2(2), dr(3), theta, rr
logical :: flag

R = radius + 1 ! radious of NP

if (size==number_NP) then
   NP = ID + 1 
else
   call random_number(rand)
   NP = int(rand*number_NP) + 1
   CALL MPI_BCAST(NP, 1, MPI_INTEGER, master, MPI_COMM_WORLD, IERR)
   NP = NP + ID
   if (NP>number_NP) NP = NP - number_NP
endif

DO
   flag = .true.
   pos = NP_position(:,NP)

   call random_number(rand)
   theta = rand * 2.0d0 * Pi

   call random_number(rand)
   rr = rand * NPmove

   x = int(rr*cos(theta))
   y = int(rr*sin(theta))
   
   pos(1) = pos(1) + x
   pos(2) = pos(2) + y
   
   ! apply PBC
   Do i = 1, 3
      if (pos(i) > dime(i)) pos(i) = pos(i) - dime(i)
      if (pos(i) < 0)       pos(i) = pos(i) + dime(i)
   enddo
   ! Check for overlap
   Do i =1,number_NP
      if (i /= NP) then
         dr = (pos - NP_position(:,i))
         dr = dr - dime * anint( dr / dime )
         dr = dr * dr
         if (sqrt(sum(dr))<2*R) flag = .false.
      endif
   enddo

if (flag) exit
enddo
NP_position(:,NP) = pos

endsubroutine move
end subroutine movenp
